home *** CD-ROM | disk | FTP | other *** search
/ Aminet 2 / Aminet AMIGA CDROM (1994)(Walnut Creek)[Feb 1994][W.O. 44790-1].iso / Aminet / util / misc / Fudgit233.lha / Source / src / medfit.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-12-14  |  2.5 KB  |  122 lines
  1. #include <math.h>
  2. #include <stdio.h>
  3. #ifndef NOUNISTD_H
  4. #include <unistd.h>
  5. #endif
  6. #include "dalloca.h"
  7. #include "head.h"
  8.  
  9. #define ERRR (-1)
  10. #ifndef NULL
  11. #define NULL 0
  12. #endif
  13.  
  14. static int ndatat=0;    /* defining declaration */
  15. static double *xt=0, *yt=0, aa=0.0, abdevt=0.0;    /* defining declaration */
  16. static double rofunc(double b);
  17. static void sort(int n, double *ra);
  18.  
  19. void Ft_medfit(double *x, double *y, int ndata, double *a, double *b, double *abdev)
  20. {
  21.     int j;
  22.     double bb,b1,b2,del,f,f1,f2,sigb,temp;
  23.     double sx=0.0,sy=0.0,sxy=0.0,sxx=0.0,chisq=0.0;
  24.  
  25.     ndatat=ndata;
  26.     xt=x;
  27.     yt=y;
  28.     for (j=1;j<=ndata;j++) {
  29.         sx += x[j];
  30.         sy += y[j];
  31.         sxy += x[j]*y[j];
  32.         sxx += x[j]*x[j];
  33.     }
  34.     del=ndata*sxx-sx*sx;
  35.     aa=(sxx*sy-sx*sxy)/del;
  36.     bb=(ndata*sxy-sx*sy)/del;
  37.     for (j=1;j<=ndata;j++)
  38.         chisq += (temp=y[j]-(aa+bb*x[j]),temp*temp);
  39.     sigb=sqrt(chisq/del);
  40.     b1=bb;
  41.     f1=rofunc(b1);
  42.     b2=bb+((f1 > 0.0) ? fabs(3.0*sigb) : -fabs(3.0*sigb));
  43.     f2=rofunc(b2);
  44.     while (f1*f2 > 0.0) {
  45.         bb=2.0*b2-b1;
  46.         b1=b2;
  47.         f1=f2;
  48.         b2=bb;
  49.         f2=rofunc(b2);
  50.     }
  51.     sigb=0.01*sigb;
  52.     while (fabs(b2-b1) > sigb) {
  53.         bb=0.5*(b1+b2);
  54.         if (bb == b1 || bb == b2) break;
  55.         f=rofunc(bb);
  56.         if (f*f1 >= 0.0) {
  57.             f1=f;
  58.             b1=bb;
  59.         } else {
  60.             f2=f;
  61.             b2=bb;
  62.         }
  63.     }
  64.     *a=aa;
  65.     *b=bb;
  66.     *abdev=abdevt/ndata;
  67. }
  68.  
  69. static double rofunc(double b)
  70. {
  71.     int j,n1,nmh,nml;
  72.     double *arr,d,sum=0.0;
  73.  
  74.     ADVECTOR(arr, 1, ndatat, "rofunc", Ft_catcher);
  75.  
  76.     n1=ndatat+1;
  77.     nml=n1/2;
  78.     nmh=n1-nml;
  79.     for (j=1;j<=ndatat;j++) arr[j]=yt[j]-b*xt[j];
  80.     sort(ndatat,arr);
  81.     aa=0.5*(arr[nml]+arr[nmh]);
  82.     abdevt=0.0;
  83.     for (j=1;j<=ndatat;j++) {
  84.         d=yt[j]-(b*xt[j]+aa);
  85.         abdevt += fabs(d);
  86.         sum += d > 0.0 ? xt[j] : -xt[j];
  87.     }
  88.     return(sum);
  89. }
  90.  
  91. static void sort(int n, double *ra)
  92. {
  93.     int l,j,ir,i;
  94.     double rra;
  95.  
  96.     l=(n >> 1)+1;
  97.     ir=n;
  98.     for (;;) {
  99.         if (l > 1)
  100.             rra=ra[--l];
  101.         else {
  102.             rra=ra[ir];
  103.             ra[ir]=ra[1];
  104.             if (--ir == 1) {
  105.                 ra[1]=rra;
  106.                 return;
  107.             }
  108.         }
  109.         i=l;
  110.         j=l << 1;
  111.         while (j <= ir) {
  112.             if (j < ir && ra[j] < ra[j+1]) ++j;
  113.             if (rra < ra[j]) {
  114.                 ra[i]=ra[j];
  115.                 j += (i=j);
  116.             }
  117.             else j=ir+1;
  118.         }
  119.         ra[i]=rra;
  120.     }
  121. }
  122.